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I. INTRODUCTION 


The nuclear reactor with delayed neutrons and tempera- 
ture-dependent feedback is a nonlinear system, whose response 
to uniform and local disturbances differs greatly from that 
Of a linear reactor. The prompt temperature feedback model 
and one average group of delayed neutrons were incorporated 
in this work. Practically all neutron dynamics analysis to- 
day deals with two-dimensional geometry, implying a symmetry 
in neutron dynamics behavior during transient [1,2]. This 
Symmetry assumption, however, could be unrealistic in the 
Safety analysis of nuclear reactors. A unique feature of 
this work is the consideration of the three-dimensional] depen- 
dence of the neutron flux. No symmetry assumptions were im- 
posed on the problem. This work investigated neutron dynamics 
under uniform and local disturbances in the core. A uniform 
initial flux throughout the interior of the reactor was im- 
posed. The finite element method (FEM) was employed to solve 
this nonlinear initial-boundary value problem. The FEM is 
quite effective in handling discontinuous forcing functions 
thereby making it particularly suited for examining the ef- 


fects of localized perturbations and space-dependent feedbacks. 
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IT. THE NUCLEAR REACTOR WITH TEMPERATURE 
DEPENDENT FEEDBACK AND DELAYED NEUTRONS 


pee HE PHYSICAL SYSTEM 

The system under consideration is a fast reactor of cy- 
lindrical geometry that is composed of two different regions. 
The reactor core or region I is cylinderical in shape and is 
fueled by U-235. Region II or the reflector region complete- 
ly surrounds the core and is composed of U-238. Both regions 
were assumed to be homogeneous. Table [I lists the physical 
properties and geometry for each region. A schematic of the 
fast reactor geomenry 1S shown in figure 1. 

Temperature and delayed neutron effects were taken into 
account. Also, a one-velocity or one-group model was assumed, 
thereby making the velocity independent of spatial or tem- 
poral effects. The delayed neutrons were considered only in 
region I since region II was assumed to be a non-multiplying 
medium. The temperature effects were also assumed to be only 
in the core region. 


In general form, the one-velocity neutron diffusion equa- 


een 15 | 3 | 

aN(r,t) = 2nN _ | + Yr 

= vDV7N - Z.vN + S(r,t) (1) 
where Sal ea Rod) 


D = neutron diffusion coefficient 
Z. = macroscopic absorption cross-section 


N(r,t)dV = number of neutrons in a volume element dV at 
a point r at time t 


Ly 





TABLE I 


Physical Constants 


DEFINITION 
radius of Region I 
total reactor radius 


height of Region I 
total reactor height 


neutron speed 


core neutron diffusion coefficient 


reflector neutron diffusion 
coefficient 


core neutron absorption cross 
section 


reflector neutron absorption 
cross section 


number of neutrons per fission 
critical fission cross section 


delayed neutron fraction; 
dollar reactivity 


fission energy 

modified convection heat 
transfer coefficient 
temperature coefficient 


abundance-weighted mean decay 
constant 


2 


VALUE 


60 cm 


90 cm 


u, 


ies 
1. 


ie 


01401 cn 


0040 cn 


005736 cn 


160 cm 
220 cm 


gx10/ cm/sec 


913 3am 
2008 em 


i 


ib 


~o4 


al 


00642 


652x101? 


cal/fission 


Oe 


0632 9 


cal/(cm*sec °C) 


-~0.004/°C 


O.42350 aaa 





a. top view 


x 
a poe ee 
2a 
I - core b. front view 
II - reflector 
Figure 1. Schematic of cylindricalmpeacron: 
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vDV*NdV = number of neutrons diffusing into dV per 
unit time at time t 


uVN dV = number of neutrons absorbed in dV per 
unit time at time t 


S(r,t)dV = number of neutrons produced in dV per 
unit time at time t 
The neutron number density is related to the neutron flux by 


the expression [4] 

o(r5t) = VNiCr,t) Oe) 
where (r,t) is the flux at time t . The neutron diffu- 
Sion equation in terms of the flux is depicted by 


i a = Dy oe - Zo eS; ate (ay) 


B. PROMPT AND DELAYED NEUTRONS 

The source or production term in equation (3) is composed 
of the contributions of the prompt and delayed neutrons. The 
majority of the fission neutrons are prompt neutrons that 
appear almost instantaneously (within ior! second) on fission. 
Assuming a fast neutron non-leakage probability of unity, the 


prompt neutron source is described by 


S (r,t) = (1-8)K (r,t) (r,t) (4) 
where 
K(r,t) = infinite medium multiplication factor 
8 = total fraction of delayed neutrons 
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The aiceiion of delayed neutrons is very small (note 
that 8 = 0.00642 from Table I). However, they have a very 
Significant effect on the reactivity because their mean life- 
times are long. Without delayed neutrons, reactor control 
would not be eciee These delayed neutrons are born in 
the decay by neutron emission of nuclei produced following 
the B-decay of certain fission fragments. For example, the 


ae plus a 


B-decay of the fission fragment Bro! leads to Kr 
neutron. Nuclei such as tate whose production in fission 
eventually leads to the emission of a delayed neutron are 
called delayed neutron precursors [4]. 

There are six main groups of delayed neutrons. Each 


group is classified according to its decay constant. fhe de- 


layed neutron source term is portrayed by 


where 


> 
it 


decay constant of the va group 


density of the gen precursor 


>) 
~_— 
“31 
ct 
a 
i 


Assuming that the fission fragments do not migrate appre- 
ciable distances and assuming a non-circulating fuel reactor 
[3], the precursor density is delineated by 

ae) 
aaah = 6.K 2.9 - ce (6) 


where B is the fraction of delayed neutrons of the jen 
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group. The solution to the precursor equation is in terms 


of a time integral expressed by 


c 
C.(#,t) = ita f eri lt-t ke CR tye(Rtydt! (7) 


0 
Inserting equation (7) in equation (5) yields the delayed 


neutron production term as 


For convenience the six main groups of delayed neutrons 
were considered as one group. This was accomplished by us- 


ing the abundance-weighted mean decay constant defined by 


5 ] 6 


This iS a reasonable approximation since many of the impor- 
tant phenomena of nuclear reactor dynamics can be character- 
ized satisfactorily by combining all the emitters or precur- 
sors into one, two, or, at most, three effective groups [3]. 
Replacing ds with the abundance-weighted mean decay con- 


Stant showed the delayed neutron source term to be 
S4(f,t) = pic, ff ett dy (a tole t)dt! (10) 


C. DOPPLER EFFECT AND TEMPERATURE FEEDBACK MODEL 
The Doppler effect in fast reactors is due to the tempera- 
ture broadening of many closely spaced high-energy resonances 


in both the fission and parasitic-absorption cross-sections. 
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These resonances basically mean that as the temperature in- 
creases, the number of neutrons that are absorbed increases 
and, similarly, the number of fission events increases. 
These nonproductive and productive processes compete in a 
complicated manner, and the net effect may be either an in- 
crease or decrease in reactivity [3]. For most fast reac- 
tors, the effect is negative, and it will be so assumed in 
this work. 

The reactivity change with respect to the fuel tempera- 


ture change is modeled by [2] 


where a,b, and c are parameters determined from experimen- 
Gal or neutronic calculations, m is an integer, and 7 1s 
the current fuel temperature. For most fast reactor cores 
(ceramic-fueled cores), TdKk /dT is very close to being con- 
stant over a wide range of temperatures. Therefore, it is 
assumed that a and c are zero and the Doppler coefficient 


is expressed by 


The following initial conditions were used: 


Eo Ko (12a) 
T(7,0°) = T, (12b) 
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The reactivity model is then expressed by 


k(t) = Ko + b en (aL) (13) 
0 
where 
KY = multiplication factor at steady-state 
vy = number of neutrons produced per fission 
Ti = original fuel temperatures 
b = Doppler constant 


Hne definition of xe is 


5 
¥ 
Kees (14) 
a 
where 
* 
Z = critical fission cross-section 
ag = macroscopic absorption cross-section of the 


Core 


The rise in fuel temperature is depicted by 


This is also described by the integral [5] 


6(7,t) =f f(t-t') w(r,t)dt! (16) 


where f(t-t') is the feedback kernel and is dependent upon 


the type of temperature feedback model used. w(r,t) is the 
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rise in neutron flux above the steady state and is expressed 


by 
w(r,t) = (r,t) - @ (r) ET) 
where d5(r) is the neutron flux at steady state. 


There are three types of temperature feedback models that 
can be considered here. The first is known as “Newton's 
feedback" which determines the reactor temperature by Newton's 
law of cooling. The second temperature feedback model is 
called the adiabatic feedback model. This represents the 
temperature for the loss of coolant case. The third model is 
the prompt temperature feedback model in which the fuel tem- 
perature follows the behavior of the neutron flux without de- 
lay [5,6]. The prompt feedback model was employed in this 
analysis. The feedback kernel for this temperature feedback 


model is 


eee) I= men cite) (us) 


<([~ 


where K is an energy production operator with units of °C 
per unit flux and y , a dimensionless quantity, is related 
to the mean time for heat transfer to coolant. 

Inserting equation (18) in equation (16) and performing 


the integration produced a rise in temperature of 


a(Ft) = 5 (Ft) (19) 
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From equation (19) the temperature of the fuel is 
= aK - : 
T(#,t) = EB y(rt) + 1h(F) (20) 


The ratio of the current temperature to the steady-state tem- 


perature is therefore 


ae ae, al (20) 


In this work Tie was considered to be constant throughout 
the core. Incorporating equation (21) into equation (13) 


gave the reactivity model of 


kK (r,t) = K° + b an (—\ wlF.t) +1] (22) 
co fo) sales 


Deweee tee EQUATIONS 

Before establishing the field equations it is desirable 
to express equation (3) in terms of the rise in flux. Using 
equation (17) in equation (3) and grouping terms yielded the 


following diffusion equation: 


te 
tL op = 2 } / -\(t-t') f 
Toe 7 LDV7y-Zy +(1-B)K2.¥ + BAZ, f e Kvdt'] 
O 
co ) 
2 x, -A(t-t' i 
+ [DV79-£9,+(1-B)K,2,0,+ ate, f e Kodt'] 


ay) 
The second bracketed term in equation (23) is identically 


equal to zero since it is the steady state portion of the 
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diffusion equation. The rise in neutron flux above its 
Steady state value is therefore expressed, for the core, by 
1 ow _ 
Vi ra = DV*yp - yy + hee Kt aU 


E t =i +! 
+ BAL, J ej ee Kv dt’ (24) 


) 


amid tor the reflector region, by 


ay _ 2 
a DV-y - ZY (25) 


<|- 


Inserting the reactivity model into equation (24) yielded 


Tow . - S = e) 
ee DV-w a (1-8) Kou 


E 
K wy = Swen!) ; 

+(1-8)br_ {2Qn[— + 1]} wv + BAK? e wdt 

a ls a 

iE y t 
+ BiEb tl eg A(t-t Ton (Se #09) deen (26) 
0) 
0 


fommeene reflector, the last four terms of edtlation | 2e mane 
zero. The effects of the temperature on the delayed neutrons 
were neglected in this work. The field equations can now be 
expressed, for the core, as 

ay 2 z s 

= vDV7y + [vz v(1-8) vi ely 
+ [-(1-8)vz_b] [en(Se + 1)]y 

a var 


+ [-BivEev] fe Be ihe) = (27) 


Zu 





and, for the reflector, as 
P) 
st - vDV*y + vE.y = 0 (28) 


The non-linear terms of the core field equation will be 
Yinearized accordingly. In more compact form, equations (27) 


and (28) became 


t a i 
st Sere aint C2) + cafin (St 6 SLi aeget )ydt']=0 
9) 0 
(29) 


and 


ye clv2y + c2y = 0 (30) 


where the meanings of the coefficients cl, ¢€2, etc. j are 
obvious for the core and reflector. 

Equations (29) and (30) were subjected to the following 
Gonaic1ons: 

boundary condition: 


Wp(7,,t) = 0 (31) 


where rp are coordinates of points on the outer surface of 
the reflector and the subscript R_ refers to the reflector 
eomcinuity of flux: 


Wel(r,t) = p (ry>t) (32) 


where ry are coordinates of points on the core-reflector 


interface and the subscript c refers to the core. 


Ze 





lil... FINE Ce BEMeh a 


A. INTRODUCTION 

The application of the finite elements of nonlinear con- 
tinua has been mostly in the field of solid mechanics. 

Prior work [1] using the finite elements of nonlinear con- 
tinua on a nonlinear reactor dynamics problem has been suc- 
cessful. The FEM in this work was utilized to reduce a non- 
linear partial differential equation of the nuclear reactor 
to a system of nonlinear ordinary differential equations in 
time. The time integration was accomplished by using a com- 
puter program for the numerical solution of stiff differen- 
tial equations developed by Franke [7]. In order to minimize 
computer storage requirements, an optimum compacting scheme 
(OCS), described by Ref. 8, was adopted. 

The finite element models of operator equations are gen- 
erally classified into three categories: a) the variational 
finite element models such as the Ritz method, b) the weigh- 
ted residuals method such as the method of Galerkin, and c) 
the direct finite element models which are not based on 
functional minimization. From experience in structural 
mechanics, the most effective method for generating accept- 
able finite element models of nonlinear equations is the 
Galerkin method [1]. This work adopted the method of Galer- 
kin in a finite element approximation over the spatial 


domain of the field equations. 
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B. THE METHOD OF GALERKIN 

The Galerkin method is a special case of the method of 
weighted residuals. It involves a rational choice of 
weighting function that is consistent with the type of fi- 
nite element approximation considered. Indeed, the weight- 
ing functions chosen are the basis or shape functions em- 
ptoyed in the finite element approximation. There are two 
favorable characteristics of the Galerkin method which 
makes it attractive. The first attribute is its amenability 
to integration by parts. This supplied the freedom of using 
a lower order finite element than might be otherwise possible. 
The second favorable characteristic of the Galerkin method 
is that the symmetric operators in the field equations trans- 
form into symmetric matrix operators. Both these attributes 
are attractive for computational purposes. 


Consider the initial-boundary-value problem 


Jt = Lv =) tr ieee) (ese) 


where af contains the nonlinear operators. According to 
the spatial finite element discretization, the solution of 
equation (33) is in the form of an union of an N-term approx- 


imation given by 


w(w,t) = o(m,t) = UW oi(t) Gy(F) , 0=1,2,...,8 (34) 


where N is the number of system degrees of freedom (ji.e., 
number of coordinates), and Gj(r) are the system or global 


basis functions which span the space of the approximate 
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oa = 


Seiution wir,t) [1]. The Einstein summation is usted. The 
global basic functions are “pyramid functions", each of 
which has a prescribed functional description over a sub- 
domain of the system and is zero elsewhere [9]. The unknown 
coordinate functions y3(t) are the time-dependent magni- 
tudes of the approximated flux w(r,t) and/or its deriva- 
tives at discrete nodal points [10]. 


The residual function, R(r,t) , is defined such that it 


~~ aS 


Hseidgentical to zero when w(r,t) is equal to the exact solu- 
tion. The residual function is expressed by 


a(7,t) = 2 -{v - f (35) 


The Galerkin orthogonality condition (using the basis func- 
tions as weight functions), when applied to the residual func- 


tion, requires that 
G,R(r,t)dVol = 0 , 1212 een icy, 


Vol 


From the field equations, the residual function for the core 


R(r,t)= oy - clv2y + c2y + cA Lin (st + el 


O 
t = 
+ c5[ / eta ) pdt' ] bey | 
6) 
and for the reflector 
Meret) = ay cly*y + c2u (38) 
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Using equation (34) and applying Galerkin's orthogonality 


condition produced a core equation of 


J ¢,16,0,(2)-996,(c1-v), + Gy(c2-y), 


Vol 
K 
tLen(or” Gyby (3) | Gy (c4eyp) | 
t -_ 
tf e~A(t-t (c5+y), dt']G 4} = 0 (39) 
0 
where RP ae 


Cy 
iT 
7 

vi 
NO 

ve 

vi 
= 


The reflector has a similar equation without the nonlinearity 


and is expressed by 


gee - VG (cl-¥), 1° Gy(c2-) jy} avols= 20 (40) 
Vol 


C. THE ELEMENT 

A three-dimensional quadratic isoparametric element was 
employed in this work. The parent element is a triangular 
prism or solid wedge with straight sides. The element shape 
functions are expressed in terms of area coordinates in the 


plane of the triangle and by an isoparametric coordinate 
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along the prism axis. Figure 2 shows the parent element. 
This element was chosen because of the ease with which it 
Mipoeene Cylindrical structure when itis transformed into 

a curved element. This type of element has been used before 
as filler elements [11]. 

The area coordinates are defined by area ratios. Con- 
Sider the triangle shown in figure 3. An arbitrary point P 
within the triangle defines three subareas designated by A,» 
A» 3 and Ax. The ratio of each of the subareas to the total 


area 1S known as an area coordinate. In equation form, the 


area coordinates are 


L, = A,/A (41a) 
Lo = An/A (41b) 
L. = A./A (ale) 


where A 1s total area of the triangle. Ly > Los and L. 
are the natural coordinates for a triangle. The requirement 
that the sum of the subareas be equal to the total area is 


obviously satisfied by the identity 
Le tk. + (a ere) 


In the plane of the triangle, only two of the area coordi- 
nates are independent. 

The isoparametric coordinates are best visualized by con- 
sidering the rectangular prism shown in figure 4. Isopara- 
metric coordinates are normalized coordinates such that 
their values on the faces of the rectangle are + 1. The nz 


axes are in general not orthogonal. They are orthogonal 
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ma 
gure 2. Quadratic triangular prism parent element 
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Figume 3. Definition of area coordinates 
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Figure 4. Isoparametric coordinates 
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only in the special case of a rectangular prism element 
[12]. The element basis functions use the ¢ coordinate 


in the prism axis and the L Los and L 3 coordinates in 


] 3 
the plane of the triangle. 

The parent element, as shown in figure 2, has 15 local 
nodal points around its periphery. As such, there are 15 


basis functions which are given below [11]: 


Sonmenr nodes: (nodes 1, 3, 5, 10, 12, 14) 


N, = 5 L,(2L,-1) (14) = FL, (1-2?) 
Na = y by (2Ly-1)(145) = 5 Ly (1-6?) 
Ne = gba (2L4-1)(14e) - $L3(1-2?) 
Naot by (2by-1) (i-s) = 5 by (8?) 
N15 P L,(2L,-1) (1-3) - 5 L,(1-¢*) 
Nyg= gp bg(2bg-1) (1-6) - 5 b3(1-0) 


Midside nodes of rectangles: (nodes 7, 8, 9) 


Pes eaGi-c°) 


7 
= _72 
= _72 
Ng = L401 e* ) 


Midside nodes of triangles: (nodes 2, 4, 6, 11, 13, 15) 
ig = ll riley 


2 ] 

Ny = ~ ll 
No = 2L4L,(1+2) 
Ni45 2L,L,(1- -o) 
N13 2L5L,(1- -¢) 
N15* 2L4L,(1- -z) 
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The coordinates of each local- node, in terms of 1 L L 


ge eee 
and ct are listed in table II. These element basis functions 
< N > define the geometry of the element. Note that they 


satisfy the relationship 


1 , at node 7 


0 , at Other nodes (43,) 
On the element level, the variation of the unknown func- 
tion w is approximated by 


O = Net eeee (44) 


where - au = row vector of element shape functions 


{w}° = column vector of time-dependent nodal 
values of w° 


To satisfy continuity requirements, the shape functions 
< N'> have to be such that the continuity of the unknown 
function wp is preserved in the parent coordinates [11]. 
The shape functions < N > which characterize the element 
geometry and the shape functions < N'> which describe the 
unknown function do not necessarily have to be the same. 
There is no requirement that the nodal values be associated 
with the same nodes which were used to define the element 
geometry, though in practice it 1s often the case. Consider 
for example, the illustrations in figure 5, If the nodes 


defining the element geometry and the nodes defining the 
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TABLE [1 


Coordinates of Local Nodal Points 


Local Node 


a 


a 


ae 


a 


ou 


ww 


10 


11 


ni 


2b 


WZ 
is 
14 
1S) 


a 


rit 


wt 


Ny 
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(a) Isoparametric 


(b) Super-parametric 


(c) Sub-parametric 
e - geometry nodes 


© -syvariable funct ionmaodes 


Figure 5. Element classification 





unknown function are identical, the element is known as an 
isoparametric element. This means that the shape functions 
describing the geometry and the shape functions describing 


the unknown function w are equal, or 
< NS =< ie > (45) 


If there are more nodes defining the geometry than nodes de- 
fining the variable function, the element is called a super- 
parametric element. Using more nodal points to define the 
unknown function than to describe the geometry of the element 
results in a Ssubparametric element [11]. This work utilized 


the isoparametric element classification. 


WeepevaviSlON OF THE SYSTEM INTO ELEMENTS 

In three-dimensional space, the division of the sys- 
tem into discrete elements is difficult to visualize. I[t is 
virtually impossible to show every nodal point of the system 
in one schematic. To present a clear view of the discretized 
domain, a "layer" approach was adopted. The first finite 
element grid or mesh employed here consisted of 128 elements. 
Under this grid (mesh I), the reactor was divided into four 
layers as shown in figure 6. Each layer was composed of 32 
elements. The first and fourth layers were each 30 cm in 
height and each contained entirely reflector elements. The 
second and third layers were each 80 cm in height and together 
they encompassed the entire core plus the remaining reflector 
elements. Each layer, in turn, was partitioned into three 


horizontal (xy) planes. The top plane included all the global 
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Layer 3 II 





I - core II ~ reflector 


Figure 6. Layers of mesh I 
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or system nodes corresponding to local or element nodes 1 
through 6. The middle plane contained all the system nodes 
corresponding to the local nodes 7, 8, and 9. System nodes 
corresponding to element nodal points 10 through 15 comprised 
the bottom plane. To fix ideas, the three nodal planes of 
the first layer of mesh I are shown in figures 7, 8, and 9. 

In this work there was only one curved side which was in 
the plane of the triangle, as shown in figure 10. The first, 
Seventh, and tenth element nodes were each arbitrarily assigned 
to have aS an opposite side the curved side of the triangle. 
The remaining local nodes in each of the respective planes 
of the element were then numbered consecutively in the counter- 
clockwise direction. 

At this point it iS appropriate to define connectivity. 

The connectivity of an element is a row vector array that 
relates the local nodes to system nodal points. The connec- 
tivity lists the system nodal points that are "connected" to 
form the element domain in the sequence of local nodal number- 
ing. Thus, the connectivity matrix for mesh I is a matrix 

of size 128 x 15. To illustrate, the connectivity of element 


number 106 1s 
eee, 210, 193, 194, 195,e211, 266, 176, 177, 283, W525 105 ieee 


A second finite element mesh (mesh II) consisting of 192 
elements was developed. It contained the same number of ele- 
ments per layer as mesh I. However, mesh II has six layers. 


The second and third layers of mesh I were each divided in 


oy 





264 233-234 


240 
256 


= 
ee) 
a] 
feed 
© 
Aas 


257 241 


256 242 


243 





Number - element number 


Figure 7. Top nodal plane of the first layer of mesh I, 
2 =r 0 cm 
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268 


Zoo 





Number - element number 


Figure 8. Middle nodal plane of the first layer of 
mesh I, z = 95 cm 
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312 281 282 


305 


304 290 


303 aol 





296) 967 206 


Number - element number 


Figure 9. Bottom nodal plane of the first layer of 
mesh I, z = 80 cm 
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Local nodal numbering of curved = 


elements 
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half, thus forming the two additional layers of mesh II. 
The connectivity matrix and the coordinates of each global 
node for mesh I are given in Appendix A. Appendix B lists 
the connectivity matrix and nodal coordinates of mesh II. 
The reader can construct the different layers and elements 
without much difficulty by using the nodal coordinates and 
the connectivity matrix. A mesh generator was not utilized 


in this work. 


E. COORDINATE TRANSFORMATION 

The use of a quadratic or higher order element permits 
the transformation or mapping of the straight-sided parent 
element into an element with curved sides. Distorted or 
curved elements provide a better fit to curved domains than 
linear elements, and thus a smaller number of elements is 
required to represent the structure adequately. 

The transformation from cartesian coordinates to curvi- 
linear coordinates can be accomplished by employing a one-to- 


one correspondence defined by [11] 


X g L, 

y |} =f Un} orf {re} (46) 
é 

Z ‘a a 


The element shape functions are utilized to achieve this 


transformation via the relation 


(47) 


<< 
1} 
—- 4 + M4 

ra 
<< 
w 

af 
i 
—s 
Ld 

nN 
= 
a) 


N 
—e [J 
me 
aw | 
N 
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where n° is the number of element nodes, and the element 
Shape functions N. are ame terms of local coordinaces. eaten 
set of local coordinates corresponds to only Onley set Oh gedies 
tesian coordinates. 

In performing the transformation, the comaare eae re- 
quirement must be met. The transformation into the new, 
curved elements should leave no gaps between adjacent elements. 
If two adjacent elements are generated from parents in which 
the element shape functions satisfy continuity requirements, 
then the curved elements will be contiguous [11]. For the 
1soparametric element, uniqueness of coordinates ensures com- 
patibility. Continuity is assured when adjacent elements 
are given the same sets of coordinates at common nodes. 

Since the element shape functions are in terms of local 
coordinates, an element of volume, dxdydz, must be transformed 
into an element of volume expressed in local coordinates. 

This is achieved through the use of the Jacobian matrix de- 
fined below. Using the chain rule, the relationship between 


E,n,t and a corresponding set of cartesian coordinates 


ne 
aN i FP az aN 
al Tees OG an 
NiN [ax ay) az i ia) 
dn on om =e y 
Nid fax ay az] | 2hs 
Ya OC goon ce dZ 
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iresJsacobian matrix [J] is defined as 


a— ? OF ° BE 

- | 2x ay az 
ax ay az 
GG ew wag 7 aie 

From equation (47), the Jacobian matrix becomes 
ON. ON. ON. 
Feo ae 8 pe Be 
; 1 ; 1 ; 1 
ON. ON aN. 
i - 3 Jie ieee Sn ae? e on oa 

j j j 
ON oN. oN. 
ie ee 


The determinant of the Jacobian is used to transform the vol- 
ume of element in cartesian coordinates to local coordinates. 


For a volume of element [11], 
dxdydz = det [J] d&dndz ole) 


Note that the determinant of the Jacobian is a variable for 
elements of curved geometry. Only in the case of straight- 
Sided elements is the determinant of the Jacobian a constant. 
In the plane of the triangle, the area coordinates (Ly; 
L,,L,) number one more than the cartesian coordinates (x,y). 
ess L is defined as a dependent variable. This establishes 


the origin of the &n coordinate system at corner point 3, as 
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illustrated in figure 11. Recall that the En axes need 


not be orthogonal. As such 


f= Ly (52a) 
nN = L. (52b) 


Using equation (42), 


L, =]-€&- 7 ( 52 e) 


oe se ( dSia) 


Me eg ey as (53b) 


Using equations (52a), (52b), and (52c) in equations (53a) 
and (53b) gives 


_ ee 
3 3L, ” aL, (54a) 
JN. oN. ON. 

mm" 3, ~ 3, (545) 
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Figure 1]. n& coordinates in a triangle 
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Now the Jacobian matrix can be evaluated using the shape func- 


tions eo N (Ly sL,,L,,¢) . Inserting equations (54a) and 
(54b) in equation (50) produces 


lt ar an =i oN. oN. 
~ ae 5 a WZ 
E aL, ol. 1 , ol, dL. 1 3L, dL. 1 

ey . * an — an oN. oN. 
9 “ie iz » Seas a y 9 x ( Hz 
lee 3 OL, aL. j aL, ols j aL, aL. j 

oN. oN. oN. 
eae Sry; 55) 


(Sib) 


To summarize, suppose it is required to transform the 


integral 


| = J F (L, L5.b452) dxdydz (56) 
Vol 


to an integral entirely in terms of local coordinates. The 


determinant of equation (55) is then utilized to give 


I fief F'(L, Lo sb,.¢)det{d(L, .L,,L,,c) Jd, diode 
-l 9 0 Gra 
Equation (56) is now in a form suitable for numerical inte- 
gration. 
The development of coordinate transformation to this 
point has been under the general assumption that all the 


sides of the element are curved. In this work, the only 
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eieved™side is in the plane of the triangle. The sides of 
the element along the prism axis are straight. As such, a 
modified Jacobian matrix can be employed, thereby reducing 
the number of calculations to be performed. This modified 


Jacobian is the 2x2 matrix defined by 


N aN. 2N aN 
( Hie pet Le Ne: et Ls 
‘ Ly aL, i aL, aL 1 
_ IN. ON. aN aN. 
[J ] = Ti) oe abe a (58) 
| OL Bla a) Se OG cL Yj 


Along the prism axis or g direction, it can be shown that 


dz = dz (59) 


nm|=> 


where h = height of the element. The volume relationship 


is then given by 


dxdydz = 3 det{J ] dL,dLjz (60) 


The integral of equation (56) then assumes the form of 


1 oT de-L 
] 
_h | * 
r= hf ff Fu ayatgotddetto"(L, sLytgoc)]dL, sdbysde (61) 
-l 0 Oo 


Equation (61) is the basis of numerical integration applied 


in this work. 


F. CONSTRUCTION OF ELEMENT MATRICES 
The system matrix operators can be constructed through 


the use of the global basis functions G, or through the 
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application of the element shape functions. Although the 
global basis functions were used to demonstrate the method 
of Galerkin, the construction of the system matrices in this 
work was achieved through element considerations. The solu- 
tion of the unknown variable wW within the element domain 


was approximated by 
oS 
= 2 Ne Wet) at T= Vee, «noe (62) 


where n° is the number of element nodal points, N. are 


the element shape functions, and we (t) are the time-depen- 
dent nodal magnitudes of we . The element contribution to 
the system matrix operators is defined by Galerkin's ortho- 


gonality condition expressed by 
/ N, R© dVol = 0 , j=1,2,...,n° (63) 


where R° is defined by replacing w with he jn equations 
(37) and (38), and the integration is over the element volume. 
Using equation (62), the element contribution is portrayed 


for the core by 


C fy, N. dol ]ty.~ it N.V2N.dVol]{ (cl + “) HE f NN. dVolJ{(c2-y") ,} 
vol Vol ° 
+ Lf NN,2m(1 + A EN, vy o(t))dVol} {(c4-v®) .} 
Vol ok 
+ [ f N,N;dvo1] { f ett!) (c5-y(t') )dt'} = 0 (64) 
Vol 
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Mere 1,J,k = 1,2,...,6 , and cl, c2, @tc., "are Constaits 
at node i. The bracketed expressions represent Square 
Matrices of size 15x15, and the braced expressions represent 
column vectors of size 15xl. For the reflector, the last 
two terms of equation (64) are zero. Before any operation 
can be performed on equation (64), the nonlinear terms (last 
two terms) must be "linearized" and the term with the 72 
Operator must be integrated by parts. 15 

The nonlinear feedback term, 2&n[l1 + ee Nv (t)] ; 
was linearized by using predicted values ls 
of the unknown function at time t. These predicted values 
come from the time integration scheme by Franke [7]. 
Basically, the integration scheme utilizes a predictor-cor- 
rector method which predicts values of the unknown function 
at the next time by using the derivatives of the function. 
Adopting the predicted values a enabled integration over 
Space. The term in equation (64) involving the nonlinear 


feedback can be written as 


[ f NN. &n(1 + TAN yy Novy’ + Scat Ny c¥ye.))dVOT] {(c4-v").} 


The last tern of equation (64) describes the delayed 

neutron contribution. In general, 
x tn x y tha r 

ety f ait va (t')dt' = Re Oe! fe tn-1 f et vec(t')dt'] 

0 6 

a) nN = \ 
te ty f ast wan (t!)dt' (65) 
t 

n- | 


where n is number of time steps, ce 1s the currentetime: 
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and th] is the previous time. To approximate the integrals 


in equation (65), the predicted values vs were employed 
in a simple trapezoidal rule. The trapezoidal rule was be- 
lieved to be sufficient since very small time steps were uti- 


lized in this work. For example, at time ts ; 


(SUM, ) = (previous SUM. )e7*" ae (eo "w.P(t,)+ W 


where 
H = current time step = to - ty 
(SUM ) i -Aty fe RSs Slt!) dt! 
‘ e W 
O 
- cy = 
(previous SUM, ) = Baa if a batt! )dt 
O 


The previous sum is formed by accumulating the trapezoidal 


integration of each time step. In general then, for time th 


(SUM, ) = (previous SUMa) 1c Jey.) (66) 


1 n-1- 


The delayed neutron term in equation (64) can be expressed 


as 


Cf N,N.dVo1] {(c5+SUM) ;} 
Vol 


Su 





In order to bring equation (64) to final form, the VY? 
operator was integrated by parts. Since the flux at the sur- 
face of the reactor is zero, integration by parts yields 


ON, ON; ON; AN, aN, aN, 
Jo N,V?Nidxdydz = -f (~t+od4 Sy yt ae ae axdydz (67) 


Vol Vol 


In vector notation, equation (67) is expressed as 


ON. 
—) 
OX 
: aN, ON oN. aN 
\ N.V?Ndxdydz = -f “say BE? Lae Pp dxdydz (68) 
Vol Vol 
N. 
* 
dZ 
Using the chain rule produces 
wnt a s 
OX OX OX dL) 3 
ea = all ui 0 aul Mi) 
ay vo oye aL. ~ aly 
ON. ON. 
1 dt 1 


Letting the 3x3 matrix above be [B'], the V? term becomes 


oN. ON. 
f (ac > ate) 
NV Ni dxdydz = 1 3 
Vol aN, ON, ON, ON, ON, ; ON, ON, 
a f Pa. e 3, - a, yaa ] [B J ‘aL, - OL,” dxdydz 
Vol 
ON. 
_ 1a (70) 
oa 


oy 





where rai]! is the transpose of [B'] . By applying the 
chain rule in equation (47) and using equation (59), it can 


be shown that for this work 





© 
© 
>|r 


rB']’ can be derived from equation (71). 
Note from equation (64) that there are three basic ele- 
ment matrices which are defined as follows after applying 


equation (60): 





h ] ] T-L, a 
(6.J-afff N.Nidet[J ]dL,dLpdz (72) 
-l 00 
, Jd i-b, aN, aN, aN, aN aN. “4 
(GG. . ] = > fff “oe, = at) * "St, = BL, aE eee ] 
-loo 
aN. ON. 
(apt - ae 
1 3 
IN. ON. r 
cB" ] es aC) det(J JdL,dL,dz 
ON 
oe (73) 
epee : : : . 
(G66, |] : aes NN on (1+ NY +...4Ny ob, —) )detld ]dLy dL ode 
O 
-loo0 


re) 
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Element matrices [G55] and [GG] are independent of 
time. However, [GGG ,; ] is time dependent due to the utili- 
zation of the predicted values v. which changes with time. 


In terms of these three basic element matrices, equation (64) 


is 
(651005) + (665, ]CcT-v"),} + [65 ,1¢(c2-v") 5) 
+ (666, .] {(c4-y") .} Se iSO SI) = 0 (75) 


The last two terms of equation (75) are zero for the reflector. 


G. CONSTRUCTION OF THE SYSTEM MATRICES 

The 15x15 coefficient element matrices were calculated 
secoraing to equations (72), (73), and (74); and “the results 
were collected element by element into the corresponding sys- 
tem coefficient matrices. The system coefficient matrix 
[BIGG] is developed from [65,1 » LBIGGG] from L6G, 51 ; 
and [BIGH] from [GGG ., J . £BIGG] and [BIGGG] are inde- 
pendent of time and can be constructed once and for all from 
geometry considerations. [BIGH] is dependent on both geom- 
etry and time due to the time dependence of the predicted 
flux utilized in the feedback term. Thus, [BIGH] is recal- 
culated at each time increment. 

Non-zero contributions to a global nodal point I come 
only from adjacent elements sharing that same nodal point I. 
Thus, the system matrices are sparse and banded. The pro- 


cess of assembling contributions from element matrices 
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Bemeires the identification of a local nodal point (i1=192), 
.,15) with a global nodal point (1=1,2,...NUMNP, where 
NUMNP is the total number of system nodes). This correspond- 
ence between element and global nodes is accomplished via 
meVemconmectiVity matrix. 
The formal treatment of the field equations in terms of 


the system coefficient matrices is described by the equation 


eae Gi {v,} + [BIGGG] {(cl-v),} + [BIGG] {(c2-v),} 


+ [BIGH] {(c4-v) 7} + [BIGG] {(c5+SUM) 7} ant) (76) 


where the system matrices are NUMNP x NUMNP and the column 
vectors are of length NUMNP x 1. However, the direct appli- 
cation of equation (76) requires a large amount of computer 
Storage. Jo take advantage of the sparsity of the system 
matrices, an optimum compacting scheme described by Ref. 8 
was employed. 

The concept behind OCS is simply to store only the non- 
zero terms of a coefficient matrix. OCS requires two integer 
arrays, say JB and NAME, and a vector of non-zero coeffi- 
cients of the square system matrix. For purposes of illus- 
tration, the square system matrix is called B, and the vector 


Mt integer 


of non-zero coefficients of B is called BB. The i 
entry in the NUMNP x 1 JB vector is the number qa: Vis 


number is defined by 


i=] 


Sc ae a Py i=1,2,...,NUMNP (775 
Jet 
where Pj is the number of terms in the jen equation. In 
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other words, Pp; 1s the number of nodes that the a node 


"sees". JB is therefore a pointer vector of length NUMNP+1 
whose 7 term locates the initial position in the BB vector 


ormecne contributing coefficients to the aa equation. The 


NUMNP 
NAME vector of length Mx1, where M = x Pas consists of 
1=] 
NUMNP successive vector blocks of variable length Pas 


Tige,..-sNUMNP. The p. integer numbers in the tel. block 


of NAME list the Ps contributors to the ql: equation. The 
Mx] BB vector contains the real non-zero coefficients of the 
NUMNPxNUMNP B matrix, arranged in the same contiguous block 
arrangement as the NAME vector. The gee term in the pee 
block or BB(JB(I) + J-1) is B(I,K), where K = NAME(JB(I)+Jd-1). 
To illustrate, consider the grid shown in figure 12. The two 


array vectors and the coefficient vector matrix of non-zero 


terms are: 


fem =) <5, 10, 13. 18, 25. 30, 33. 38, 42> 
NAME = <1,2,4,5/2,3,1,6,5| --- |9,8,6,5> 
BB = <By 1 2By 5 »By 4 2B 5 | Boo »Bo3 »Boq Bog Bos| ---| Bag »Bgg »Bgg Bag” 


In this illustration, NUMNP = 9 and M = 41 [8]. 

In this work, a judicious method of numbering the system 
nodal points was adopted to further reduce computer storage 
requirements. Since the surface boundary nodes of the reac- 
tor represent zero neutron fluxes, the contributions of these 
nodes to interior or non-zero nodes can be discarded. Thus, 
only the interior nodal points need to be considered. These 


non-zero nodes were numbered first in the finite element mesh 
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2 8 


Number - element number 


Figure 12. Sample grid used for illustrating OCS 
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used so that in the OCS, the number of non-zero nodes (NNZ) 
replaces NUMNP. 

The vectors of non-zero coefficients will be designated 
BIGG, BIGGG and BIGH since the square coefficient matrices 
described in equation (76) were not utilized. To illustrate 
the application of OCS in the system, the following sample 


program is given: 


Boee500 J=JBB, JE 
LL = NAME(J) 
DY(I1) = DY(I) + BIGG(J)*W(LL) + BIGGG(J)*cl1(LL)* 
W(LL) + BIGG(J)*c2(LL)*W(LL) + BIGH(J)*c4(LL)* 
W(LL) + BIGG(J)*c5(LL)*SUM(LL) (78) 
500 CONTINUE 
450 CONTINUE 


where 
LL = nodal point "seen" by node I 
JE-JBB = total number of nodal points that node I "sees" 
DY(I1) = summation of all contributions to node I! and 
should sum to zero as stated by the field equa- 
tion 
The assumption of a homogeneous reflector was relaxed in the 
application of equation (78). Interface nodes were assigned 


properties of the core. Therefore, if LL is a reflector 
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node not on the core reflector interface, the last two terms 


of equation (78) are nonexistent. 


59 





IV. NUMERICAL INTEGRATION 


A. LINE AND AREA INTEGRATION 

The solution of the element matrix equations was achieved 
through numerical integration since an exact closed form solu- 
tion cannot be established. The volume integration was accom- 
plished by using a line integration in the g-direction and an 
area integration in the plane of the triangle. The line inte- 


gral is described by the Gaussian quadrature formula [12] 
] 


n 
Jf scree = 2, H,f(a,) (79) 
oy k=] 


where nis the number of Gauss integration points 


BLS 
UN 


weighting coefficients 


the function f(z) evaluated at Gauss point a, 


mame III lists ays Hy eure (cpm | El 


The area integration was achieved by the equation 


1 a = 
2 m, m, om 
{f F(L,,L,,l,)dL dL, * s Wei 5 ote) (80) 
0 0 


m= | 


where m is the number of area integration points and w. 
are the weights. The numerical integration points for the 


area integration are given in Table IV which was extracted 
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PxAeee Tit 


Abscissae and Weight Coefficients 
of the Gaussian Quadrature Formula 





] n 
ao mtede = Spey Hs ayy) 
| k=] 
+a H 
n= 2 
0:57735 02691 89626 1-00000 00000 00000 
. n=3 
0:77459 66692 41483 0:55555. 55555. 35556 
0:00000 00000 00000 088888 88888 88889 
n=4 
0:-86113 63115 94053 0-34785 48451 37454 
0:33998 10435 84856 0:65214 51548 62546 
nz 5 
0:90617 98459 38664 0:23692 68850 56189 
0:53846 93101 05683 0:47862 86704 99366 
0:00000 00000 00000 0:56888 88888 88889 
n= 6 
093246 95142 03152 0°17132 44923 79170 
0:66120 93864 66265 0°36076 15730 48139 
0:23861 91860 83197 0:46791 39345 72691 
a) 
0:94910 79123 42759 0:12948 49661 68870 
0:74153 11855 99394 0:27970 53914 89277 
0-40584 51513 77397 0°:38183 00505 05119 
0-:00000 00000 00000 0:41795 91836 73469. 
n= 8 
0:96028 98564 97536 0:10122 85362 90376 
0:79666 64774 13627 0:22238 10344 53374 
0:52553 24099 16329 9°331370 66458 77887 
0:18343 46424 95650 0°36268 37833 78362 
n=z=9 
0:96816 02395 07626 0:08127 43883 61574 
0:83603 11073 26636 0:18064 81606 ©4857 
0:61337 14327 00590 0°26061 06964 v2935 
0-32425 34234 03809 0:31234 70770 40002 
0:00000 00000 00000 0:33023 93550 9!260 
n= 10 
0:97390 65285 17172 0:06667 13443 08688 
0:86506 33666 88985 0:14945 13491 50581 
0:67940 95682 99024 6-21908 63625 !5982 
0:43339 53941 29247 0:26926 67193 09996 
0:14887 43389 81631 0:29552 42247 14753 
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Numerical Formulas for Triangles 


TABLE IV 





Order Fig. Error Points Tnangular Werghts 
Co-ordinates 2W, 
Linear AY R = 0{h’) a 4.4.4 
a 4, 4,0 { 
Quadratic R = 0(h3) b 0,4.4 4 
c= 4, 0,4 § 
a f.4.4 -h 
Cubic AN R = QA‘) b Wak ra] 
=—_ 4 co Bt 
d at 
a 4444 35 
b 44.0 
¢ 0.4.4 \ 46 
Cubic os R = O(A*) d  4,0.4 
g ae e 1.0.0 
} 0.1.0 } a 
g 0,0. 1 
a t4.4 0:225 
'\ b Xie Py 
c B,. x1. By 0-13239415 
Quintic A> R = Mh") d By. By. 1 
RN e 2Z2- 82, 
<< % f Bae In Br 0-12593918 
g 
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with 

2, = 0-05971587 
B, = 047014206 
B, = 010128651 





PrOtmeanre?. il... The volume integration of the function 


F(L,.L,.L,,8) using numerical integration is therefore 


eal 
Pf Sf sey rgetgecraargae : 
0 O 


n m 
Dy ace! ea (81) 
oA it i a 
where Wy Hy In the manner of equation (81), the element 
matrices become 
h n m m m m 2K = 
[G..J = 5 uk ay we N 5 (Ly Lo sL3 .C )N.det[d I (82) 
h n m m m m _k * 
(6G, 5] ai 2 Wr » We F (Ly ees 33 ) detf{J ] (83) 
= m= | 
h n m m m m _k ‘ 
(ace, ;] = > ye Wie 2 Me gl (Lol GES ce) eecer| Os limeces 
k=] m= | 3 


where F “and T are easily derived from equations (73) and 
(74), respectively. Note that N, and det[J ] are also 
evaluated at each integration point. They are not shown as 
such merely for the sake of convenience in writing the equa- 


ETONS . 


B. NUMBER OF INTEGRATION POINTS 

It is difficult to estimate the number of integration 
points required for good accuracy due to the complexity of 
the functions involved. The basic rule that the best number 
of integration points is found by trial and experience was 


adopted. For the cre element matrix, the numbers involved 
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can be approximated. This was done by letting the determi- 
nant of the Jacobian equal to twice the area of the curved 
triangle (checking det[J*] at each integration point showed 
that this assumption was not too unreasonable). With det[J*] 
outside the integration process, [G54] can be solved in 

closed form by integrating out ¢ from -1] to +1 and then apply- 


ing the closed form equation [12] 


| 
Mm, mM. m my:m,im I 


3 " 3 
L, L, L. d(Area) = 2A Cm +m54m42) 1 
Area 
(85) 
where My, M5, m, are positive integer exponents and A_ is 


the area of the triangle. 

Five test points within the 15x15 element matrix were 
selected, as listed in Table V. The values obtained from 
the application of equation (85) are also given in Table V. 
Three sets of area integration points, each with a different 
number of ¢ Gauss points, were used. These three sets of 
area integration points, given in Table IV, are: 

1. e¢Wbic order (4 points) 

2. cubic order (7 points) 

3. quintic order (7 points) 

Using each of the three integration points above with differ- 
ent ct Gauss points in equation (82) produced the results 
obtained in Table V. From these results, the quintic order 
area integration points were selected for the element matrix 
re with three c¢ Gauss points in the prism axes. This 


set of integration points was also used for (GGG, ]. 
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TABLE V. 


Selection of Integration Points for LG. 5] 


Area 
@apoilnts oints G(1,] G22 G(1,9 G(6,9 G(9,9 


ie cubic 1415 5278 -2756 5072 8536 
(4 pts) 

5 : 1817 5278 -3170 5072 10240 

7 : 1817 5278 -3170 5072 10240 

3 cubic 3248 8247 -3114 4960 10243 
nese) 

5 , 3249 8247 -3114 4960 10240 

7 ; 3249 8247 -3114 4960 10240 

3 quintic 2464 6600 -3144 5020 10240 
(7 pts) 

5 , 2464 6600 -3144 5020 10240 

7 y 2464 6600 -3144 5020 10240 

Approximated values 2500 6700 -3142 5026 10053 
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The numbers for the [GG..] element matrix could not 
be approximated. The best that could be done was to obtain 
an idea of the order of magnitude of this element matrix. 
Towards this end, a linear triangular element was assumed. 
Using linear approximation, the order of magnitude was found 
to be about 1g Using the three different area integration 
points mentioned above with varying ¢ Gauss points in equa- 
tion (83) yielded the results given in Table VI. Further 
enecks on the [GG 5. ] element matrix showed that the cubic 
order with four area integration points yielded the desired 
order of magnitude. Thus, the fourth order cubic with five 
> Gauss points was employed for the (6G... ] element matrix. 
A note should be mentioned here in regards to the vast dif- 
femences in results obtained for (GG. .] using different 
area integration points. Most likely, it was due to the 
fea) and rai]! matrices which required the inversion of 
x, ; =* , etc. Or perhaps it was caused by the nature of 
the hybrid element used in this work. In any case, further 


investigation is warranted in this area. 
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C 


3 


Oints 


TABLE VI 


Selection of Integration Points for (GG. | 


Area 
oints 


CUDI1C 
(4 pts) 


quintic 
(7 pts) 


GG(1,9) 


18.4 


22.0 
22.0 
-1.03x10 


-1.03x10 
-].03x10 
-52.6 


-52.6 
-52.6 


4 


4 
b 


GGC3.3 
uel 


L776 


va 
2.11x107 
2.11x107 


211x107 


Apel re 


1.27x10" 
1.27x10° 


6/7 


-1.05x10 


-6.64x10 
-6.64x10 


GG(14,9 
-165 


- 186 
- 186 


1 OSealg? 


7 
1.05x107 


4 


4 
4 


GG(8,8 
16529 


1959 


195.9 
8.42x10/ 
8.42x10/ 


8.41x10" 


6.84x10° 


6 84x10" 


> 


6.84x10 


GG Sede 
106.1 


106.1 
106.1 
2056 


2056 
2056 
1.51x10 


TeSPxct0 
-1.51x10 


5 


5 
5 





V. TEST PROBEEMS WWD. RESULTS 


The reactor was subjected to uniform and local perturba- 


tions in the form of a ramp input described by 
ee ee Fae (86) 


Mere & is the change in Le per unit time and a is 
me critical fission cross-section. oa must first be ob- 
tained before the perturbations can be applied. This was ac- 
complished by trial and error until a stationary solution 


* 
was reached. For meshI, 2 was found to be 0.005/360 


¢ 
per cm. No attempt was made to find uf” for mesh II due 
to time limitations. As such, the test problems outlined be- 
low were applied to mesh I. Future work is planned to apply 
the test problems on mesh II. 

The following perturbations were applied: 


a) Uniform perturbation of 10 dollar of reactivity per 


second: 
Ze(r,t) aan at , in the core 
where a = 0.005893/cm-sec 


b) Local perturbation at the core center of 100 dollar 


of reactivity per second: 


* - = 
P +atd(r) , in the core 


68 





where Py 1s) (050,0) and o = 0.015407/em-<ec 
jm focal, off-center perturbation: 


* =~ : ° 
Ze(r,t) pele. a.té(r,) > i = 1,2338% in themcere 


where 
ry = (0,60,40) 
a, = 0.015407/cm-sec = 100 dollar per second 
Go = 0.008123/cm-sec = 50 dollar per second 
a, = 0.005894/cm-sec = 10 dollar per second 


Three test points, (0,0,0), (60,0,0), and (-60,0,80), 
were selected to trace the neutron time history. For cases 
a) and b), the neutron flux was plotted at each test point 
during transience. This is shown in figures 13 and 14. 

Case c) involved three ramp inputs and were conducted for 
both the linear and nonlinear reactor equations. The linear 
and nonlinear responses were compared at each test point for 
each ramp input and are illustrated in figures 15 through 23. 
The radial and axial flux distributions at time t = 0.0123 
Second were plotted for the steady state and the 100 dollar 
perturbation case. These are embodied in figures 24 and 25. 
Finally, a neutron flux early time history between mesh I 
and mesh II was plotted to check the effect of using a finer 
element mesh. The result is portrayed in figure 26. 

Figures 13 and 14 revealed a clear space dependence of 
the neutron flux during transience as expected. Figures 15 


through 23 demonstrated the effects of temperature feedback 
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neutron flux 
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10 
1011 
101? 
10° 
10> 
A 
B 
6 
Erdure 13. 


ne ioe 10°* 


time (sec) 


neutron flux at test point (0,0,0) 
neutron fiUx ae test pomt oo. ond 
neutron flux at test point (-60,0,80) 


Neutron flux transient history at three 
test points with a uniform perturbation 
of 10 dollar of reactivity per second. 
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neutron flux 





time (sec) 


A - neutron flux at test point (0,0,0) 
B - neutron flux at tesuapoint (60,0,0) 
C - neutron flux at test.point (-60,0.,80) 


Figure 14. Neutron flux transient history at various 
test points for a local central perturbation 
of 100 dollar/sec of reactivity. 
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Figure 15. Linear and nonlinear fluxes at (0,0,0) due 
to a local perturbation of 100 dol lar/sec een 
reactivity at (0,60,40). 
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A - nonlinear neutron response 
B - linear neutron response 


Figure 16. Linear and nonlinear fluxes at (60,0,0) due 
to a local perturbation of 100 dollar/sec at 
(Om ou. 40) . 
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A - nonlinear neutron response 
B - linear neutron response 


Linear and nonlinear fluxes at (-60,0,80) 
due to a local perturbation of 100 dollar/sec 
at (0,60,40). 
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neutron flux 





Ene (sec) 


A - nonlinear neutron response 
B - linear neutron response 


Figure 18. Linear and nonlinear fluxes at (0,0,0) due 
to a 50 dollar/sec local perturbation at 
Oo. 40.) 
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Figure 19. Linear and nonlinear fluxes at (60,0,0) 
due to a 50 dollar/sec local perturbation 
Gree Oo 40. 
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Figure 20. Linear and nonlinear fluxes at (-60,0,80) 
due to a 50 dollar/sec local perturbation 
ae 0,00, 40). 
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Figure 21. Linear and nonlinear fluxes at (0,0,0) due 
to a 10 dollar/sec local perturbation at 
Heo 60,40) . 
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neutron flux 





times seer 


A - nonlinear neutron response 
B - linear neutron response 


Figure 22. Linear and nonlinear fluxes at (60,0,0) due 
to a 10 dollar/sec local perturbation at 
(0,60,40). 
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neutron flux. 


Firaure 23. 





10> 10° Mee len 


Time seer 
A - nonlinear neutron response 
B - linear neutron response 


Linear and nonlinear fluxes at (-60,0,80) 
due to a 10 dollar/sec local perturbation 
ac Om 0ee4.0 )).. 
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60 30 0 228, ~60 


x=0, z=0, yy 


@ - point of local perturbation 
A - neutron flux at steady-state 
B - neutron flux under local perturbation 


figure 24. Radial flux distribution for thessteaqy stage 
and 100 dollar/sec local perturbation. 
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pagure 25. Axial flux” distribution for thessteady =stapeacag 
100 dollar/sec local perturbation. 
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A - mesh I 
B - mesh II 


Figure 26. Early time history of the neutron flux at core 
center using mesh I and mesh II. 
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and delayed neutrons on the flux. To obtain an idea of the 
difference in magnitude between the linear and nonlinear 

mlx tor the 100 dollar local perturbation, at time t = 107° 
second, the linear case predicted a neutron flux at the core 


center of 2.024 x 10!¢ neutrons per cmé per sec. The non- 


linear reactor predicted 1.85 x one neutrons per em® per 
sec. Although the numbers are not sufficiently accurate due 
to the crudeness of the finite element mesh, the significant 
difference in the order of magnitude between the linear and 
nonlinear neutron flux leads to the belief that the three- 
dimensional finite element model utilized is predicting the 
correct trend of neutron flux behavior. 

The radial flux distribution for the steady state shown 
in figure 24 portrays the expected flux distribution. The 
radial flux distributions for the local perturbation case 
Show a skew distribution at the point of perturbation. The 
axial flux distribution at steady state 1S very much symme- 
tric about the center. Again, the expected skew at the 
point of perturbation is there, as shown in figure 25. 
Figure 26 gives an indication that the oa) for finer 
meshes is higher. However, curve B of figure 26 should be 
extended to longer times to verify this hypothesis. Curve B, 
as plotted in figure 26, used four hours of computer time 


employing the H-compiler of the IBM 360/67. 





VI. CONCLUSIONS AND RECOMMENDATIONS 


The finite element mesh employed here is crude, and, 
thus, results that were obtained should be considered as posi- 
tive indicators rather than numerically conclusive facts. The 
trend of neutron flux behavior is the major thrust of this 
work. From the results obtained, it is concluded that the 
expected patterns of neutron behavior as predicted by the 
three-dimensional finite element model used do occur. These 
patterns were best demonstrated by the differences between 
the linear and nonlinear flux responses and by the spatial 
flux distribution at steady state and during local perturba- 
tion. The three-dimensional quadratic finite element model 
utilized should produce better results by resorting to a 
finer mesh. The draw-back to a finer mesh, of course, 1s 
the significant increase in computer time and storage re- 
quirements. Once accurate results are obtained through 
finer element meshes, comparisons between three- and two- 
dimensional models can be attempted. 

A mesh generator was not developed for this problem. It 
is recommended that this be done to ease the transition from 
One mesh to another and to minimize human error. In addi- 
tion, a similar calculation using a three-dimensional linear 
element should be performed to corroborate the results ob- 
tained here. It should be particularly noted that this type 


of problem is highly sensitive to the fission cross-section. 
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* 
iaeune searcn for Le > it was necessary 


to adjust 2 to 


‘9 


the sixth decimal place. There is no exact method of deriv- 


* 
ing Ue due to the highly nonlinear aspect of the problem; 


therefore, trial and error must be used. 

Gaussian quadrature used to determine the 
Should be investigated. The cause of the 
values obtained by using different number 
points should be established. Was it due 
points, the coordinate transformation, or 
functions? This is an important question 


numerical results could be based. 
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Finally, the 
element matrices 
differences in 

of integration 

to the integration 
the element shape 


upon which future 





APPENDIX A 


MESH I CONNECTIVITY AND COORDINATES 


follows 


1S aS 


The connectivity matrix for the 128-element mesh 
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The coordinates of the 128-element mesh are as follows 
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APPENDIX B 


MESH II CONNECTIVITY AND COORDINATES 
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